### This script will compare empirical measurements of microclimate to three climate predictions
### The three predictions are BIOCLIM, AWAPCLIM, and microCLIM
### Should also plot daily microclimate for daily microCLIM predictions

### Load libraries

library('SDMTools')

### Set working directories

out.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/Plots/AWAPvsEMP/'

### Read in Raw Data file with preds

raw.data = read.csv('/home1/99/jc152199/brt/FINALOPTIMALMODELS/Model_Data_Plus_Preds.csv',header=T)

### Read in min model data

#load('/home1/99/jc152199/brt/FINALOPTIMALMODELS/OptimalMinModel.Rdata')

### Summarize raw data to monthly values by site
### Start by appending a column for month

raw.data$month = as.numeric(substr(raw.data$date,5,6))
raw.data$year = as.numeric(substr(raw.data$date,1,4))

### Next aggregate micro_max and micro_min

### Use a for loop to aggregate by site

allsitemonths.out = NULL # Blank object to bind data onto

for (site in c(unique(as.character(raw.data$site))))

	{
	
	# Subset raw data to a single site
	
	t.data = raw.data[which(raw.data$site==site),]# & is.na(raw.data$micro_max)==F & is.na(raw.data$micro_min)==F),]
	t.data = na.omit(t.data)
	
		for (yearx in c(unique(t.data$year)))
		
			{
			
			# Subset t.data to a single year
			
			tt.data = t.data[which(t.data$year==yearx),]
	
				for (monthx in c(unique(tt.data$month)))
		
					{
					
					# Subset tt.data to a month year
					
					ttt.data = tt.data[which(tt.data$month==monthx & is.na(tt.data$micro_min)==F & is.na(tt.data$micro_max)==F),]
					
					# Check to make sure the month has data in it before summarizing
					
					if(nrow(ttt.data)!=0)
					
					{
					
					# Calculate monthly meanmin and meanmax
					# Preds 1 are from the original optimal BRT model
					# Preds 2 are from the BRT model with all vars except rainfall and realized radiation and longitude
					
					ttt.out = data.frame(site=site,year=yearx,month=monthx,obs_meanmax=mean(ttt.data$micro_max),obs_meanmin=mean(ttt.data$micro_min),preds_meanmin=mean(ttt.data$min_preds),preds_meanmax=mean(ttt.data$max_preds),AWAP_meanmin=mean(ttt.data$AWAPmin),AWAP_meanmax=mean(ttt.data$AWAPmax),n.dayspermonth=nrow(ttt.data))
					
					# Assemble data to write out, only if the month has more than 10 days of data
					
					if(ttt.out$n.dayspermonth>=10) {allsitemonths.out = rbind(allsitemonths.out,ttt.out)}
					
					}
			
					cat(site,'-',yearx,'-',monthx,'\n')
					
					}
			
			}
			
	}
	
# Write out summary data frame to all site year months	
	
write.csv(allsitemonths.out,file=paste(out.dir,'Monthly_MeanMax_&_MeanMin_From_Empirical_Data.csv',sep=''),row.names=F)
			
### Perform summary stats on the dataframe allsitemonths.out

monthlysummary = aggregate(allsitemonths.out[,c(4:9)],by=list(site=allsitemonths.out$site,month=allsitemonths.out$month),FUN=mean)

### Use a loop to calculate the monthlymean temp from monthly meanmax and meanmin in the dataframe monthlysummary

monthlymean=NULL # Blank object to bind data onto

for (i in c(1:nrow(monthlysummary)))

	{
	
	t.data = monthlysummary[i,]
	
	tt.data = data.frame(site=t.data$site,month=t.data$month,obs_meantemp=mean(c(t.data[,3],t.data[,4])),preds_meantemp=mean(c(t.data[,5],t.data[,6])),AWAP_meantemp=mean(c(t.data[,7],t.data[,8])))
	
	monthlymean = rbind(monthlymean,tt.data)
	
	}
	
# Close loop

t.bc01 = aggregate(monthlymean[,c(3:5)],by=list(site=monthlymean$site), FUN=mean) # Mean Annual Temperature	
		
t.bc05 = aggregate(monthlysummary[,c(3,6,8)],by=list(site=monthlysummary$site), FUN=max) # Max Temp of Warmest Period
			
t.bc06 = aggregate(monthlysummary[,c(4,5,7)],by=list(site=monthlysummary$site), FUN=min) # Min Temp of Coolest Period

monthlymean.K = data.frame(site=monthlymean$site,month=monthlymean$month,obs_meantemp.K=monthlymean$obs_meantemp+273.15,preds_meantemp.K=monthlymean$preds_meantemp+273.15,AWAP_meantemp.K=monthlymean$AWAP_meantemp+273.15)

t.sdofthemean.K = aggregate(monthlymean.K[,c(3:5)],by=list(site=monthlymean.K$site), FUN=sd) # SD of monthly mean temps in Kelvin

t.bc01.K = aggregate(monthlymean.K[,c(3:5)],by=list(site=monthlymean.K$site), FUN=mean) # Calculate monthly mean in Kelvin

t.bc04 = data.frame(site=t.bc01.K$site,Obs_CofV=t.sdofthemean.K[,2]/t.bc01.K[,2],Preds_CofV=t.sdofthemean.K[,3]/t.bc01.K[,3],AWAP_CofV=t.sdofthemean.K[,4]/t.bc01.K[,4]) # CofV for monthly mean temps

t.bc04[,2:4]=t.bc04[,2:4]*100 # Express as a percentile
			
final.out = data.frame(site=t.bc01$site,obs_bc01=t.bc01[,2],obs_bc04=t.bc04[,2],obs_bc05=t.bc05[,2],obs_bc06=t.bc06[,2],preds_bc01=t.bc01[,3],preds_bc04=t.bc04[,3],preds_bc05=t.bc05[,3],preds_bc06=t.bc06[,3],AWAP_bc01=t.bc01[,4],AWAP_bc04=t.bc04[,4],AWAP_bc05=t.bc05[,4],AWAP_bc06=t.bc06[,4]) # Assemble data into a dataframe to write out

write.csv(final.out,paste(out.dir,'BC_Values_From_Empirical_&_AWAP_&_Predicted_Data.csv',sep=''),row.names=F)
			
#####################################################################################################	
### Now create some plots comparing AWAP and microCLIM preds to Empirical
### May need to subset by sites later

#for.subset = aggregate(allsitemonths.out$n.dayspermonth,by=list(site=allsitemonths.out$site),FUN=sum)

#complete.sites = for.subset[which(for.subset$x>365),]

#sub.final.out = final.out[which(final.out$site %in% complete.sites$site),]

#sub.final.out = final.out

#sub.ASCII.out = ASCII.out[which(ASCII.out$site %in% sub.final.out$site),]

### Small dataframe of predtypes

#preds.df = data.frame(num=c(0,4),preds=c(1,2))

outdata = NULL

for (ii in c(1,4,5,6)) # Loop through all surfaces, determining the 99% CI's of the slope and intercept for the relationship between empirical and predicted data

	{
	
	t.final.out = data.frame(site=final.out$site,final.out[,which(substr(names(final.out),nchar(names(final.out)),nchar(names(final.out)))==ii)]) # Subset data by surface type
	
	AWAP_lm = lm(t.final.out[,2]~t.final.out[,4]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	mc_lm = lm(t.final.out[,2]~t.final.out[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	
	
	t.AWAP.sigdata = data.frame(pred='AWAP',mod=mod, surface=ii, min.slope = summary(AWAP_lm)$coefficients[2,1]-(2.58*summary(AWAP_lm)$coefficients[2,2]), slope.pred = summary(AWAP_lm)$coefficients[2,1], max.slope = summary(AWAP_lm)$coefficients[2,1]+(2.58*summary(AWAP_lm)$coefficients[2,2]), min.int = summary(AWAP_lm)$coefficients[1,1]-(2.58*summary(AWAP_lm)$coefficients[1,2]), int.pred = summary(AWAP_lm)$coefficients[1,1], max.int = summary(AWAP_lm)$coefficients[1,1]+(2.58*summary(AWAP_lm)$coefficients[1,2]))
	t.mc.sigdata = data.frame(pred='mc',mod=mod, surface=ii, min.slope = summary(mc_lm)$coefficients[2,1]-(2.58*summary(mc_lm)$coefficients[2,2]), slope.pred = summary(mc_lm)$coefficients[2,1], max.slope = summary(mc_lm)$coefficients[2,1]+(2.58*summary(mc_lm)$coefficients[2,2]), min.int = summary(mc_lm)$coefficients[1,1]-(2.58*summary(mc_lm)$coefficients[1,2]), int.pred = summary(mc_lm)$coefficients[1,1], max.int = summary(mc_lm)$coefficients[1,1]+(2.58*summary(mc_lm)$coefficients[1,2]))

	outdata = rbind(outdata,t.AWAP.sigdata,t.mc.sigdata)
	
	}
	
# Close

# Calculate significance of slope, to be a non-signficant difference, the expect value (1) must fall within the 99% CI's

outdata$slope.sig = NA
outdata$slope.sig[which(outdata$min.slope>1)]='sig'
outdata$slope.sig[which(outdata$max.slope<1)]='sig'
outdata$slope.sig[which(outdata$min.slope<1 & outdata$max.slope>1)]='nonsig'

# Calculate significance of intercept, to be a non-signficant difference, the expect value (1) must fall within the 99% CI's

outdata$int.sig = NA
outdata$int.sig[which(outdata$min.int>0)] = 'sig'
outdata$int.sig[which(outdata$max.int<0)] = 'sig'
outdata$int.sig[which(outdata$min.int<0 & outdata$max.int>0)] = 'nonsig'

	

for (ii in c(1,4,5,6)) # Loop through 'BC' summaries of the different surfaces

	{
	
	t.final.out = data.frame(site=final.out$site,final.out[,which(substr(names(final.out),nchar(names(final.out)),nchar(names(final.out)))==ii)]) # Subset data by surface type
	
	AWAP_lm = lm(t.final.out[,2]~t.final.out[,4]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	mc_lm = lm(t.final.out[,2]~t.final.out[,3]) # Perform a linear model using the Empirical Data as y-values and the preds (AWAP or microCLIM) as x-values
	
	#newdata1 = data.frame(obs=t.final.out[,4]) # Bind x-data for plots into a dataframe
	#newdata2 = data.frame(obs=t.final.out[,3])
	
	#AWAP.rawpoly = predict(AWAP_lm, newdata1, interval="confidence"); AWAP.rawpoly = cbind(newdata1,AWAP.rawpoly) # Calculate 95% CI's for all x values in the AWAP linear model
	#mc.rawpoly = predict(mc_lm, newdata2, interval="confidence") ; mc.rawpoly = cbind(newdata2,mc.rawpoly)# Calculate 95% CI's for all x values in the microCLIM linear model
	
	#AWAP.polybottom = AWAP.rawpoly[c(order(newdata1)),c(1,3)]
	#names(AWAP.polybottom) = c('x','y')
	#AWAP.polytop = AWAP.rawpoly[c(rev(order(newdata1))),c(1,4)]
	#names(AWAP.polytop) = c('x','y')
	#AWAP.polyout = rbind(AWAP.polybottom,AWAP.polytop) # Calculate a polygon that represents the 95% CI's of the AWAP linear model
	
	#mc.polybottom = mc.rawpoly[c(order(newdata2)),c(1,3)]
	#names(mc.polybottom) = c('x','y')
	#mc.polytop = mc.rawpoly[c(rev(order(newdata2))),c(1,4)]
	#names(mc.polytop) = c('x','y')
	#mc.polyout = rbind(mc.polybottom,mc.polytop) # Calculate a polygon that represents the 95% CI's of the microCLIM linear model
	
	lims = range(c(t.final.out[,2],t.final.out[,3],t.final.out[,4]),na.rm=T) # Calculate limits for x and y axes
	
	png(paste(out.dir,'Surface0',ii,'.png',sep='')) # Open .png device driver
		
	plot(t.final.out[,3],t.final.out[,2],xlab='Predicted',ylab='Empirical',main = paste('Surface 0',ii,' Comparison',sep=''),col='blue',cex=1, type='p', xlim=lims, ylim=lims, pch=2) # Configure plot space
	
	#abline(b=1,a=0,col='black') # 1:1 line
	
	#polygon(AWAP.polyout[,1],AWAP.polyout[,2], col='#FF000025', lty=1, lwd=1, border=NA) # Plot a polygon representing the 95% CI's for the AWAP linear model
	
	points(t.final.out[,4],t.final.out[,2],col='red',pch=1) # Points for AWAP and Empirical
	
	abline(b=summary(AWAP_lm)[[4]][2],a=summary(AWAP_lm)[[4]][1],col='red',lwd=2) # Plot a line representing the relationship between AWAP and Empirical data
	
	#polygon(mc.polyout[,1],mc.polyout[,2], col='#0000FF25', lty=1, lwd=1, border=NA) # Plot a polygon representing the 95% CI's for the microCLIM linear model
	
	abline(b=summary(mc_lm)[[4]][2],a=summary(mc_lm)[[4]][1],col='blue',lwd=2) # Plot a line representing the relationship between microCLIM and Empirical data
		
	#points(t.final.out[,3],t.final.out[,2],col='blue',cex=.5,pch=2) # Points for microCLIM and Empirical
	
	legend('topleft',legend = c('Summarized Microclimate Model',paste('Adj. r^2 ',substr(summary.lm(mc_lm)[9],1,4),sep=''),paste('Slope ',round(mc_lm$coefficients[2],2),sep=''),paste('Intercept ',round(mc_lm$coefficients[1],2),sep='')),text.col=c('black','blue','blue','blue'), bty='n')
	legend('bottomright',legend = c('Summarized AWAP Model',paste('Adj. r^2 ',substr(summary.lm(AWAP_lm)[9],1,4),sep=''),paste('Slope ',round(AWAP_lm$coefficients[2],2),sep=''),paste('Intercept ',round(AWAP_lm$coefficients[1],2),sep='')),text.col=c('black','red','red','red'), bty='n')
	
	dev.off()
	
	}
	
